
#include"fft.h"

struct complex spinfac[65536];
//struct complex a[131072];
//struct complex b[131072];
struct complex tem[131072];

struct complex complexmult(struct complex a, struct complex b)
{
    return (struct complex){(a.re * b.re) - (a.im * b.im), a.re * b.im + a.im * b.re};
}
struct complex complexadd(struct complex a, struct complex b)
{
    return (struct complex){a.re + b.re, a.im + b.im};
}
struct complex complexrev(struct complex a)
{
    return (struct complex){-a.re, -a.im};
}

void countspinfac(){
    int i = 0;
    do{
        spinfac[i].re = cosf(pidiv65536 * i);
        spinfac[i].im = sinf(pidiv65536 * i);
        spinfac[32768 + i].im = spinfac[i].re;
        spinfac[32768 + i].re = -spinfac[i].im;
        i = i + 1;
    } while (i < 32768);
}

void fft131072(struct complex *a, int n)
{
    int step = 131072 / n;
    int i = 0, halfn = n / 2;
    struct complex a1, b1;
    do
    {
        a1 = complexadd(a[i], a[i + halfn]);
        b1 = complexmult(complexadd(a[i], complexrev(a[i + halfn])), spinfac[i * step]);
        a[i] = a1;
        a[i + halfn] = b1;
        i = i + 1;
    } while (i < halfn);
    if (n != 2)
    {
        fft131072(a, halfn);
        fft131072(&(a[halfn]), halfn);
    }
}

void changeorder(struct complex *a, struct complex *b)
{
    int i = 0;
    int ir = 0;
    do{
        ir = ((i & 1) << 16) + ((i & 2) << 14) + ((i & 4) << 12) + ((i & 8) << 10) + ((i & 16) << 8) + ((i & 32) << 6) + ((i & 64) << 4) + ((i & 128) << 2) + (i & 256) + ((i & 512) >> 2) + ((i & 1024) >> 4) + ((i & 2048) >> 6) + ((i & 4096) >> 8) + ((i & 8192) >> 10) + ((i & 16384) >> 12) + ((i & 32768) >> 14) + ((i & 65536) >> 16);
        tem[ir] = a[i];
        a[i] = b[ir];
        i = i + 1;
    } while (i < 131072);
    i = 0;
    do{
        b[i] = a[i];
        a[i] = tem[i];
        i = i + 1;
    } while (i < 131072);
}


void setnum(short *a,short *b){
    int i = 0;
    do{
        a[i] = sinf(200 * 3.1415926535 * 0.02 * i);
        b[i] = cosf(150 * 3.1415926535 * 0.01 * i);
        i = i + 1;
    } while (i < 131072);
}
/*
int main()
{
    printf("start\n");
    clock_t st1 = clock();
    countspinfac();
    clock_t end1 = clock();
    printf("计算旋转因子完成：%d\n", end1 - st1);
    st1 = clock();
    setnum();
    end1 = clock();
    printf("数据产生完成：%d\n", end1 - st1);
    st1 = clock();
    fft131072(a, 131072);
    fft131072(b, 131072);
    end1 = clock();
    printf("fft完成：%d\n", end1 - st1);
    st1 = clock();
    changeorder();
    end1 = clock();
    printf("数据转换完成：%d\n", end1 - st1);
    int i = 0;
    st1 = clock();
    do
    {
        a[i].im = -a[i].im;
        b[i].im = -b[i].im;
        i = i + 1;
    } while (i < 131072);
    fft131072(a, 131072);
    fft131072(b, 131072);
    changeorder();
    i = 0;
    do{
        a[i].im = -a[i].im/131072;
        b[i].im = -b[i].im/131072;
        a[i].re = a[i].re / 131072;
        b[i].re = b[i].re / 131072;
        i = i + 1;
    } while (i < 131072);
    end1 = clock();
    printf("逆变换完成：%d\n", end1 - st1);
    i = 0;
    do{
        printf("%d,a,%f,b,%f\n", i, a[i].re, b[i].re);
        i = i + 1;
    } while (i < 300);
}*/
